!> author: 左志华
!> date: 2022-07-16
!>
module test_sph_internal_force

    use testdrive, only: new_unittest, unittest_type, error_type, check
    use sph_kinds, only: rk
    use sph_region, only: region_t
    use sph_math, only: get_distance
    use sph_internal_force, only: internal_force_with_visc, internal_force_without_visc
    use twin_sph_list_type, only: twin_list_t, twin_list_t_init, twin_list_t_push
    use sph_eos, only: c0
    use sph_smoothed_kernel_function, only: skf
    use sph_env, only: nums
    implicit none
    private

    public :: collect_internal_force

contains

    subroutine collect_internal_force(testsuite)
        type(unittest_type), intent(out), allocatable :: testsuite(:)

        allocate (testsuite, source=[ &
                  new_unittest("FUNC: internal_force_with_visc", test_internal_force_with_visc), &
                  new_unittest("FUNC: internal_force_without_visc", test_internal_force_without_visc) &
                  ])

    end subroutine collect_internal_force

    subroutine test_internal_force_with_visc(error)
        type(error_type), intent(out), allocatable :: error
        type(region_t) :: region
        type(twin_list_t) :: pairs
        real(rk) :: dx(2), r, dwdx(2), w

        ! 部署 region 和 pairs，一个实粒子、一个虚粒子
        nums(1) = 2
        nums(2) = 1
        nums(4) = 3
        allocate (region%particles%loc(2, 3), source=reshape([0.0_rk, 0.0_rk, 1.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], [2, 3]))
        region%particles%itype = [-2, 2, 2]
        region%particles%rho = [1010.0_rk, 1000.0_rk, 980.0_rk]
        region%particles%mass = [1010.0_rk, 1000.0_rk, 980.0_rk]

        allocate (region%particles%acc(2, 3), source=0.0_rk)
        allocate (region%particles%vel(2, 3), source=10.0_rk)
        region%particles%vel(:, 2) = [0.0_rk, 0.0_rk]     ! 速度差导致粘性力
        allocate (region%particles%u(3), region%particles%c(3), region%particles%p(3))
        region%particles%hsml = 1.0_rk

        ! 计算光滑核函数的值
        call twin_list_t_init(pairs, 2, .true.)
        call get_distance(region%particles%loc(:, 1), region%particles%loc(:, 2), dx, r)
        skf%type = 'cubic_spline'
        call skf%setup_kernel(region%particles%hsml)
        call skf%kernel(r, dx, w, dwdx)
        call twin_list_t_push(pairs, [1, 2], w, dwdx)

        call get_distance(region%particles%loc(:, 2), region%particles%loc(:, 3), dx, r)
        call skf%kernel(r, dx, w, dwdx)
        call twin_list_t_push(pairs, [2, 3], w, dwdx)

        ! 部署 EoS 函数参数
        c0 = 0.1_rk

        call internal_force_with_visc(region%particles, pairs, region%particles%acc, 3)
        call check(error, region%particles%acc(1, 2), 0.37981716771311935E-4_rk, more="internal_force_with_visc::x")
        if (allocated(error)) return

        call check(error, region%particles%acc(2, 2), 0.73637181877366114E-4_rk, more="internal_force_with_visc::y")
        if (allocated(error)) return

        ! 没有速度差和压强差
        region%particles%acc = 0.0_rk
        region%particles%rho = 1000.0_rk
        region%particles%mass = 1000.0_rk
        region%particles%vel = 0.0_rk
        call internal_force_with_visc(region%particles, pairs, region%particles%acc, 3)
        call check(error, region%particles%acc(1, 2), 0.0_rk, more="internal_force_with_visc::x")
        if (allocated(error)) return

        call check(error, region%particles%acc(2, 2), 0.0_rk, more="internal_force_with_visc::y")

    end subroutine test_internal_force_with_visc

    subroutine test_internal_force_without_visc(error)
        type(error_type), intent(out), allocatable :: error
        type(region_t) :: region
        type(twin_list_t) :: pairs
        real(rk) :: dx(2), r, dwdx(2), w

        ! 部署 region 和 pairs，一个实粒子、一个虚粒子
        nums(1) = 2
        nums(2) = 1
        nums(4) = 3
        allocate (region%particles%loc(2, 3), source=reshape([0.0_rk, 0.0_rk, 1.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], [2, 3]))
        region%particles%itype = [-2, 2, 2]
        region%particles%rho = [1010.0_rk, 1000.0_rk, 980.0_rk]       ! diff: 10, 20
        region%particles%mass = [1010.0_rk, 1000.0_rk, 980.0_rk]

        allocate (region%particles%acc(2, 3), source=0.0_rk)
        allocate (region%particles%u(3), region%particles%c(3), region%particles%p(3))
        region%particles%hsml = 1.0_rk

        ! 计算光滑核函数的值
        call twin_list_t_init(pairs, 2, .true.)
        call get_distance(region%particles%loc(:, 1), region%particles%loc(:, 2), dx, r)
        skf%type = 'cubic_spline'
        call skf%setup_kernel(region%particles%hsml)
        call skf%kernel(r, dx, w, dwdx)
        call twin_list_t_push(pairs, [1, 2], w, dwdx)

        call get_distance(region%particles%loc(:, 2), region%particles%loc(:, 3), dx, r)
        call skf%kernel(r, dx, w, dwdx)
        call twin_list_t_push(pairs, [2, 3], w, dwdx)

        ! 部署 EoS 函数参数
        c0 = 0.1_rk

        call internal_force_without_visc(region%particles, pairs, region%particles%acc, 3)
        call check(error, region%particles%acc(1, 2), 0.34104630662549003E-4_rk, more="internal_force_without_visc::x")
        if (allocated(error)) return

        call check(error, region%particles%acc(2, 2), 0.68209261325098019E-4_rk, more="internal_force_without_visc::y")
        if (allocated(error)) return

        ! 没有速度差和压强差
        region%particles%acc = 0.0_rk
        region%particles%rho = 1000.0_rk
        region%particles%mass = 1000.0_rk
        call internal_force_without_visc(region%particles, pairs, region%particles%acc, 3)
        call check(error, region%particles%acc(1, 2), 0.0_rk, more="internal_force_without_visc::x")
        if (allocated(error)) return

        call check(error, region%particles%acc(2, 2), 0.0_rk, more="internal_force_without_visc::y")

    end subroutine test_internal_force_without_visc

end module test_sph_internal_force
